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Abstract: Erosion is an important issue in soil science and is related to many environmental problems, 
such as soil erosion and sediment transport. Establishing a simulation model suitable for soil erosion 
prediction is of great significance not only to accurately predict the process of soil separation by runoff, 
but also improve the physical model of soil erosion. In this study, we develop a graphic processing unit 
(GPU)-based numetical model that combines two-dimensional (2D) hydrodynamic and Green-Ampt (G-A) 
infiltration modelling to simulate soil erosion. A Godunov-type scheme on a uniform and structured 
square grid is then generated to solve the relevant shallow water equations (SWEs). The highlight of this 
study is the use of GPU-based acceleration technology to enable numerical models to simulate slope and 
watershed erosion in an efficient and high-resolution manner. The results show that the hydrodynamic 
model performs well in simulating soil erosion process. Soil erosion is studied by conducting calculation 
verification at the slope and basin scales. The first case involves simulating soil erosion process of a slope 
surface under indoor artificial rainfall conditions from 0 to 1000 s, and there is a good agreement between 
the simulated values and the measured values for the runoff velocity. The second case is a tiver basin 
experiment (Coquet River Basin) that involves watershed erosion. Simulations of the erosion depth 
change and erosion cumulative amount of the basin during a period of 1—40 h show an elevation 
difference of erosion at 0.5—3.0 m, especially during the period of 20—30 h. Nine cross sections in the 
basin are selected for simulation and the results reveal that the depth of erosion change value ranges from 
—0.86 to —2.79 m and the depth of deposition change value varies from 0.38 to 1.02 m. The findings 
indicate that the developed GPU-based hydrogeomorphological model can reproduce soil erosion 
processes. These results are valuable for rainfall runoff and soil erosion predictions on rilled hillslopes and 
river basins. 
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1 Introduction 


Soil erosion is the foundation of soil and water conservation, watershed management, and 
nonpoint pollution control. It is also closely related to soil degradation, ecological restoration, 
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dam construction, flood prediction, and water resources management. Soils need to be managed 
even as they respond to changes in rainfall. Storm rainfall is a key and dynamic factor in eroding 
soil, and most soil erosion is caused by rainstorm events. All of these factors have great impacts 
on both land use and soil erosion. Additionally, the potential impacts of climate change and water 
resources on the frequency and magnitude of rainfall in conjunction with subsequent changes in 
the form of land use are notable. Several studies have examined the relationship between rainfall 
response to erosion and the important factors that affect gully erosion change (Qin et al., 2018a, 
b). Individual heavy rainfall events can also cause severe soil erosion (Wang et al., 2016; Bryndal 
et al., 2017). Researchers have performed much work on the effects of individual rainfall events 
and rainfall intensity characteristics on soil erosion (Xu et al., 2013; Wang et al., 2016; Cerdà et 
al., 2017; dos Santos et al., 2017; Chen et al., 2018; Panagos and Katsoyiannis, 2019). Typically, 
severe soil erosion is caused by several heavy rainfall storms. The erosion process of the Loess 
Plateau in the arid region of Northwest China is very prominent, and some scholars have 
estimated the amount of erosion especially in the rill erosion process, which is very important to 
understand soil erosion in this area (e.g., Shen et al., 2016). 

Rainfall plays an important role in soil erosion and transport because soil particles formed by 
erosion are transported to form gullies that hinder the further use of land (Woodward, 1999; 
Borrelli et al., 2017). Furthermore, more water flow energy is converted into kinetic energy for 
separation, which is entrained and transported, thereby exacerbating the erosion rate (Wirtz et al., 
2010). Over the last decade, many indoor and outdoor experiments have been carried out to study 
soil erosion (González-Hidalgo et al., 2009; Taye et al., 2013; Wu et al., 2017, 2019; Qin et al., 
2018a, b; Wang et al., 2018; Zhao et al., 2019; An et al., 2020). Empirical models are commonly 
used in field scale research of soil erosion, such as the Universal Soil Loss Equation (USLE) 
proposed by Wischmeier and Smith (1978) and the Revised Universal Soil Loss Equation 
(RUSLE) proposed by Renard et al. (1997). Most of these models are calculated using a series of 
measured data, and the changes in these data can explain soil erosion. The parameters in the 
empirical USLE and RUSLE models consider only the length and slope but ignore the spatial 
characteristics between rills and slope shape. In terms of the separate problem of gully erosion, 
scholars have aimed to model the rainfall erosion process through physical mechanisms (An and 
Liu, 2009; Qin et al., 2018a, b). The Water Erosion Prediction Project (WEPP) model was 
developed on the basis of actual physical processes; it defines erosion as the sum of material 
transport between rills. Similar physical models have also been developed and applied, such as 
the Kinematic Runoff and Erosion Model and European Soil Erosion Model (Smith et al., 1995; 
Morgan et al., 1998). In addition, there are other factors affecting the evolution processes of gully 
width, such as the discharge and slope (Bingner et al., 2016). Wells et al. (2013) found that with 
increases in the slope and flow rate between gullies, the width and widening rate of rills also 
gradually increased. The collapse process and physical mechanism of the sidewall for a slope 
caused by erosion are very complicated and are of widespread interest (Chaplot et al., 2011; Chen 
et al., 2013; Bingner et al., 2016). Although much research on soil erosion has been carried out 
and fruitful results have been achieved, the damage caused by gully erosion in different regions is 
still difficult to predict through existing models, and the parameters need to be localized. At 
present, numerical simulation techniques can simulate the above processes, but there are still 
some difficulties in quantitative simulation of erosion depth (Abderrezzak et al., 2013; Wells et 
al., 2013; Bingner et al., 2016). Accordingly, further in situ experiments and model studies are 
needed to increase the understanding of ravine retreat (Chaplot et al., 2011). A better 
understanding of the gully expansion process is essential for the development of some water 
erosion prediction model modules (Wells et al., 2013; Bingner et al., 2016; de Lima et al., 2019) 
to protect erosive gullies. 

Several researchers have developed two-dimensional (2D) runoff and soil erosion models that 
can simulate the spatial characteristics of rill erosion, which can determine the rainfall-runoff 
process and simulate the spatial characteristics of slope runoff (Govindaraju et al., 1992; Gottardi 
and Venutelli, 1993; Howes et al., 2006; Chen et al., 2013; Aksoy et al., 2016). However, these 
models lack physical modules that can characterize the process of soil erosion. Therefore, 


KANG Yongde et al.: Two-dimensional hydrodynamic robust numerical model of soil erosion... 


one-dimensional (1D) and 2D coupled models were developed to simulate rainfall runoff and soil 
erosion processes on hillsides (Tayfur, 2007; An and Liu, 2009). Such models idealize the slope 
as a combination of some 1D aligned rills divided into 2D interrills. On this basis, 2D slope 
erosion models were established, but they cannot distinguish erosion between different rill 
configurations (Nord and Esteves, 2005). Although many results have been achieved in the past 
ten years, including 1D and 2D models, certain problems are not suitably solved; for example, 
there is a lack of information on distributed erosion that can reflect slope and trench erosion. 
Modelling methods may incorporate the quantitative prediction of soil erosion; however, the 
interaction mechanism of the spatial structure of rills and slope shape on rainfall and soil erosion 
processes is not clear. 

The spatial characteristics of rills and nonplanar slopes of distributed models still lack research 
at the landscape scale. Some researchers have found that both rainfall and concave hillslopes 
dramatically affect soil erosion processes (Qin et al., 2018a, b; Tian et al., 2017; Batista et al., 
2019). Furthermore, for both rill erosion and slope erosion, hillslopes can produce a series of 
depressions, which provide channels for surface runoff and sediment sources (Shen et al., 2016; 
Tian et al., 2017). The direction angle, local density, and spatial characteristics of width and depth 
of rills are crucial for runoff and soil erosion processes on slopes (Appels et al., 2011; Di Stefano 
et al., 2017, 2018; Wu et al., 2019). Slope shape is an important factor affecting slope flow and 
erosion (Sweeney et al., 2015). Different slope forms in nature may lead to various runoff and 
erosion characteristics, substantially, increasing the local rill shear stress and rill erosion (Gordon 
et al., 2012; Batista et al., 2019). Due to the complexity of slope erosion process (Stroosnijder, 
2005), rainfall is unevenly distributed and complex, but sufficient representativeness plays a very 
important role in accurate runoff forecasting (Beven, 2012). A failure to consider the complex 
changes in soil surface microtopography may lead to large errors in the soil erosion model 
according to Wu et al. (2019). However, in current laboratory experiments, the duration of rainfall 
is short, the scale of soil flume is small, and the soil surface cannot form obvious topographic 
characteristics, such as rills (Aksoy et al., 2016). In general, there is no distributed modelling 
method that can reflect complex gully networks and slope shapes and provide the errors of soil 
erosion prediction. To date, no consensus has been reached on the spatial structure of rills and 
slopes in terms of rainfall, runoff and soil erosion (Wu et al., 2017, 2019), hindering the further 
development of targeted soil conservation strategies that take spatial variability into account. 

In this study, we propose a new model, namely, a soil erosion model coupled with 
hydrodynamics and the Green-Ampt (G-A) model, for simulating soil erosion to improve the 
understanding of the interaction between rainfall and soil erosion. Theoretically speaking, water 
flow on a slope is a complex three-dimensional water flow problem. Therefore, it is not 
appropriate to apply the traditional shallow water equation (SWE) to describe the water flow on a 
steep slope. It should be noted that the model established in this study is not suitable for steep 
slopes or uneven terrain. The developed model uses graphic processing unit (GPU)-based 
acceleration technology when solving the numerical model, which accelerates the computational 
efficiency. With the use of simulations of soil erosion, transport and deposition on slopes and 
small watersheds, the spatial distribution characteristics of soil erosion caused by rainfall can be 
evaluated. 


2 Modelling approach and development 


Based on the spatial and temporal dynamic model coupling physical and empirical laws, we 
design a model including soil infiltration, erosion and sediment transport modules to provide 
directions for further solving the 2D sediment transport equation. 

We calculate the overland flow by using 2D SWEs, which include one continuity equation and 
two equations of motion of coordinate directions x and y. The rainfall energy is the essence of 
sediment detachment, and the runoff energy is the conveyer energy for eroded sediment. 
Raindrops can increase the solute transport process from soil into surface runoff, and separated 
substances may be transported because of the interaction of raindrops and water alone or in 
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combination. When suspended sediment is no longer supported by water, the sediment is 
deposited due to the decrease in water transport capacity. 

Assuming that the initial soil layer is dry, rainfall should be considered in the model. In this 
case, the flow depth, velocity and sediment concentration at the calculated point are converted to 
the boundary condition being set to zero during the calculation in the model. In addition, the soil 
surface may become dry again after rainfall events. As a result, numerical complexity requires 
specific program to resolve. One case involves a very shallow water inrush, in which the speed 
requires special treatment (Hou et al., 2013a, b). The corrected program in the next time step uses 
the wetting/drying front and necessitates implementing a negative water depth method flexibly 
according to the need (Hou et al., 2013a, b). The numerical model is represented as a square 
divided in a uniform grid with longitudinal and transversal slopes in the x and y directions of the 
physical domain. These directions are represented by the upper end with a closed boundary and 
the lower end with an open boundary. Based on the rectangular grid, we couple the sediment 
transport module on the basis of the hydrological hydrodynamic numerical model, and establish a 
2D water-sediment coupling model on the slope under the moving bed condition by using the 
finite volume method, which is used to simulate the unsteady flow of water and the transport of 
viscous sand in the process of slope erosion. 

This is a fully coupled model of water and sediment that considers the influences of sediment 
density and topography on flow movement. The central upwind scheme is used to solve the 
interface flux, and the linear reconstruction of the interface variable is incorporated to ensure 
second-order spatial accuracy. The model uses the unbalanced sediment transport formula to 
calculate the sediment transport process. In the entire calculation process, the model can ensure 
that the calculated water depth is nonnegative and has strong stability and robustness. The erosion 
rate is considered to be positive under detachment and negative under deposition. The soil 
detachment rate increases with increasing flow shear stress, upon which the flow shear stress that 
acts on the soil exceeds the critical shear stress, and rill detachment can be predicted. 


2.1 Hydrology model of soil infiltration 


The approximate stagnant water model was proposed based on capillary theory (Green and Ampt, 
1911) and was termed the G-A infiltration model (Chu, 1978). In view of the advantages of the 
simple form of the G-A model, it has been improved under specific conditions for ease of 
application. The governing equation is as follows: 


R tSt 


K,|1«(6,-6)S,/1,] t»t, 
where f, represents the instantaneous infiltration rate (cm/s); R represents the rainfall intensity 
(cm/min); K, represents the saturated vertical hydraulic conductivity (cm/min); 6 and @, are the 
initial water content of soil and saturated water content (cm"/cm ), respectively; Sis the humid 
front suction (cm); J, represents the cumulative infiltration after ponding (cm); t represents the 
time (min); and t, is the time to ponding (min). 

The function f, is defined as: 


Ij (1) 


t,=—. (2) 


2.2 Soil erosion sediment transport submodule 


As an effective approximation of Saint Venant's equation (Hou et al., 2020), the governing 
equation uses the 2D SWEs of Equation 3, including the vertical average. We assume that 
pressure distribution of hydrostatic Reynold's average, viscosity and turbulent pressure are not 
taken into account. The equation is expressed as follows to ensure the conservation of its format 
(Chen et al., 2013): 
EM + ue + ec S, (3) 
Ot Ox oy 
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where q represents the flow variable vector that consists of 7; F and G are the flux vectors of the 
conserved variables in the x and y directions, respectively; S represents the source vector; A is the 
flow depth (m); qx and q, are the unit-width discharge values in the two Cartesian directions 
(m^/s); h is the flow depth (m); c is the volumetric concentration in sheet flow layer 
(dimensionless); z, is the bed elevation (m); Ø denotes the velocity difference between 
sedimentary facies and fluid water (Greimann et al., 2008), which equals 1 (Guan et al., 2014); v 
is the average flow velocity in y direction (m/s); C is the volumetric concentration in flow depth 
(dimensionless); P is the sediment material porosity (dimensionless); $5 represents the bed slope 
in source terms; S, represents the friction slope in source terms; i is the rainfall intensity (m/s); f, 
is the infiltration rate through the soil surface (cm/s); y is a coefficient; ay is the angle between S 
and x (?); g 1s the gravitational acceleration (m/s *). D is the sediment deposition rate, and it can be 
calculated by D=@,C,, where c, represents the sedimentation velocity of sediment particles (m/s), 
and C,is the near-bed equilibrium concentration at the reference level; and EF is the sediment 
stripping rate, which can be calculated by E=a@,C,,., where Cae is the balanced near-bed sediment 
concentration at the reference height level (Smith and Mclean, 1977). 


2.3 Numerical scheme 


In this study, based on the finite volume method in the Godunov format, we use the 
Harten-Lax-van Leer contact (HLLC) approximate Riemann solver to obtain the coupled solution 
of the water flow and solute transport flux. In all the test cases considered in this work, two types 
of boundary conditions, i.e., an open boundary and a closed boundary, are used. The boundary 
conditions are controlled by flux calculations at the boundary (Hou et al., 2013a, b). The solute 
advection process simulation adopts the second-order explicit total variation diminishing (TVD) 
scheme with a limiter to control the numerical error, which reduces the numerical dissipation and 
numerical oscillation caused by the advection term of the transport equation; the implementation 
is detailed in the preliminary work (Horton, 1993; Hu and Cao, 2009; Liang and Marche, 2009; 
Hou et al., 2013a, b, 2015). 

Additionally, the running time on a specific computing device is often used to quantitatively 
evaluate the computational efficiency of a numerical model. The computing device configuration 
used is an Intel 17-6700 central processing unit (CPU) running at 3.40 GHz with a Windows 10 
operating system. The C++ (C++ Programming) and CUDA (Computer Unified Device 
Architecture) languages are used to program the GPU-based accelerated parallel computing 
process. In CUDA programming, the CPU is used as the host, and the GPU is used as the 
coprocessor. The CPU is responsible for processing logical transactions and serial operations, and 
the GPU is responsible for performing parallel computing tasks. First, the data are read into the 
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CPU, and the grid information, boundary conditions and calculation parameters are initialized in 
each grid. Then, the data are allocated to the corresponding GPU device space, and the data are 
copied to the GPU memory. After the calculation is completed, the results are synchronized with 
the CPU. 


3 Model verification 


The evaluation data of the model set are performed at the State Key Laboratory of Soil Erosion 
and Dryland Farming on the Loess Plateau, Northwest A&F University, Shaanxi Province, China 
(Liu et al., 2006). The experimental parameters of the pan include a slope length of 3.2 m, a width 
of 1.0 m, a depth of 0.3 m and a slope of 20°. The artificial rainfall is produced with a side-spray 
automatic rainfall simulation system. The nozzle height is approximately 16.0 m, the rainfall 
uniformity is greater than 90%, and the rainfall intensity could be controlled in the range of 
150-200 mm/h. In the experiment, the artificial filling method is used. Local loess with a 
thickness of 25 cm is evenly laid in the groove, and the median particle size is approximately 0.02 
mm. The volumetric density of the soil is 1.33 g/cm’, the initial soil water content is 0.2206 
cm’/cm’, and the saturated water content is 0.5027 cm’/cm’. In this experiment, the soil effective 
water conductivity is 1.67x10 ° m/s, and the matric suction is 0.15 m. When rainfall runoff is 
generated, a sample is taken every other time, and the beaker volume is 300 mL each time. The 
sampling time is recorded with a stopwatch. The instantaneous runoff flow and sediment 
concentration are determined by the drying method, and then, the cumulative runoff volume and 
erosion amount are approximated (Liu et al., 2006). 

Figure 1 represents the relationships of time with runoff velocity, flow discharge and 
cumulative erosion amounts. In the first 10 min, with the time increasing, the simulated value of 
the average flow velocity shows the same sharp increasing trend as the measured value. After 10 
min, the value gradually shows a steady trend. The average flow velocity value remains at 
approximately 0.23 m/s, which is higher than the measured value. The value is slightly larger and 
lasts until 50 min, and the simulated value and the measured value tend to be similar. The flow 
rate change also shows the same trend as the flow rate. As time increases, the flow rate shows a 
sharp increasing trend. With the accumulated sediment yield increases over time, the simulated 
value and the measured value show the same trend, which is basically the same. Overall, the 
average velocity and flow rate increase with increasing time, and finally, the cumulative sediment 
yield increases, which is in line with the physical mechanism of erosion and sediment yield. The 
simulated values of flow velocity, flow discharge and cumulative erosion are closer to the 
measured values, which shows that the experimental and simulated values have the same change 
trend, indicating a good agreement between them. 

Overall, the simulated results of the flow velocity, flow discharge and cumulative erosion show 
that the model can reproduce the changes in the rainfall runoff process and erosion, and the 
capability allows efforts of erosion modelling to also be extended to slopes by using the method. 

To quantitatively analyze the simulation results, this study uses three indicators, including root 
mean square error (RMSE), coefficient of determination (R^) and mean relative error (MRE) to 
evaluate the model effect. These three parameters can be defined as follows: 

Mia (M; -S 


RMSE = , |S 7 Ui 9 
i (9) 


(n -m)(s-5) 


i (10) 
(M; -MJ xus) 
PE (11) 


where M is the measured value; S is the simulated value; M is the mean measured value; S is 
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the mean simulated value; and n is the sample number. The closer the R? is to 1, the closer the 
result is to the exact solution. The closer the RMSE is to 0, the closer the result is to a perfect fit. 

As shown in Table 1, the proposed model is shown to be reliable in simulating slope soil 
erosion and can obtain good slope morphology prediction results. The comparative analysis 
shows that RMSE, Re and MRE can also reflect the accuracy of the simulation results. Therefore, 
the prediction model constructed in this study has strong advantages. 
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Fig. 1 Relationships between flow velocity and time (a), between flow discharge and time (b), and between 
cumulative erosion amount and time (c) 
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Table1 Performance evaluation indices of the model developed in this study 


Parameter RMSE R? MRE 
Flow velocity 0.152 96.5 1.395 
Flow discharge 0.059 94.1 1.208 
Cumulative erosion 0.037 90.2 1.433 


Note: RMSE, root mean squared error; R?, coefficient of determination; MRE, mean relative error. 


4 Results and discussion 


4.1 Numerical simulation example of an indoor slope 


The first case chosen to test this model is a slope erosion test under artificial rainfall conditions 
(Wang et al., 2018). The experiment is conducted in the Rain Test Hall of the State Key 
Laboratory of Ecological Water Conservancy Engineering in Northwest Arid Areas, Xi'an 
University of Technology, Shaanxi Province, China. The test adopts an adjustable slope steel 
channel, with the slope of 15° and the test scale of 6.0 mx1.5 mx0.6 m (lengthxwidthxdepth) 
(Wang et al., 2018). To simulate the water permeability characteristics of natural soil, we arrange 
seepage holes at the bottom of the sink. The test soil samples are taken from the local loess soil in 
Ansai District (36°45'N, 109°1 1'E) of Shaanxi Province, and the samples are then air-dried. Soil 
particle components are shown in Table 2. 


Table 2 Particle components of the experimental soil 


Particle percentage (%) 


Soil type >1.000 0.250-1.000 . 0.050—0.250 0.010—0.050 0.005—0.010 0.001—0.005 «0.001 
mm mm mm mm mm mm mm 

Loess 0.00 1.30 32.80 46.70 2.90 5.80 10.50 
Natural sand 0.00 62.18 24.60 11.95 1.27 0.00 0.00 


In the procedure, we cover the loess with 20 cm of natural sand, and control the average dry 
volume of the soil at 1.09 g/cm’. Clay loess is used to quantify the effects of runoff and sediment 
from the upper slope on downslope rill erosion and sediment under different sediment 
concentrations and rainfall intensities. We simulate continuous rainfall with a downward rainfall 
sprinkler system. The total test rainfall time 1s 60 min, and the designed rainfall intensity is 50 
mm/h. After rainfall erosion starts, sediment samples are collected under the water tank at an 
interval of 5 min. 

Under the action of rainfall and runoff, the erosion evolution of rills on slopes can be divided 
into the following three processes. First, under the action of rainfall, runoff gradually forms and 
flows along rills. Second, as the water flow transports erode soil particles, the uneroded rill 
sidewalls collapse under the action of a shearing force. Third, the erosion of the toe from the 
sidewall of the rill groove starts suddenly and unpredictably because of the gravity. In addition, 
while small narrow grooves form, cracks also form on the side surfaces of the narrow grooves due 
to the action of the shearing force. When gravity is greater than the cohesion and friction of the 
soil, small collapses occur to form rills. Martinez-Casasnovas et al. (2004) and Chen et al. (2013) 
found that under the continuous repetition of rainfall and currents, the overhanging layer 
continues to be destroyed by rapid concentrated flow, leading to the acceleration of collapse and 
gully erosion. 

The specific process is described as follows. As shown in Figure 2a, at £—0 s (where ¢ is the 
time), as the rainfall time increases, runoff starts and gradually forms, and rills also form on the 
slope. Over the period of #=10—50 s, with the gradual formation of runoff, the runoff velocity of 
the slope rapidly increases, but the runoff energy consumed by the soil erosion process gradually 
decreases at the bottom of the slope and is accompanied by the formation of rills. During the 
period of £—100—1000 s, with the gradual widening and scouring of the sidewalls of the slope rills, 
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erosion channels are formed, sand particles begin to be transported, which presents small rills one 
by one, and the shape of the entire slope changes. As shown in Figure 2, these three processes 
accordingly promote and influence one another within this time. Comparatively, the change in the 
flow depth between rills is small in terms of soil erosion. As shown in Figure 2c-d, the changes in 
space are easily detectable over the period of t=100—1000 s, mainly because the strength of rill 
erosion primarily depends on the rainfall intensity and rainfall duration. The area with strong 
slope erosion is in the middle part because of the strong erosion ability at this location. Although 
the slope length is the same, the erosion intensity varies greatly in different rills. For slopes with 
greater continuity and larger inclination angle, the runoff velocity and erosion capacity are larger. 
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Fig. 2 Accumulation of slope erosion at different times. (a), 0 and 10 s; (b), 20 and 50 s; (c), 100 and 200 s; (d), 
500 and 1000 s. 
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In the actual field environment, the length of the slope usually ranges from several metres to 
nearly hundreds of metres, and the slope changes are also different. The seasonal distribution of 
precipitation is highly uneven, and the intensity of rainstorms is also very different. Most 
rainstorms have short durations and uncertain locations (Wang and Jiao, 1996; Jiang, 1997). 
Accordingly, these extreme differences in topography and rainfall can also lead to different 
characteristics of slope soil erosion. Rainfall is often random and uneven, and the variation in 
intensity has different impacts on slope soil erosion. Some scholars have studied the influences of 
different rainfall patterns on slope soil erosion (Cerdà et al., 2017). Under the same rainfall 
intensity, with the passage of time, the permeability becomes increasingly lower and more stable, 
and the excess rainfall rate also increases. This is because when the rainfall peak is later, the 
runoff erosion capacity at the peak time is higher (Liu et al., 2006). 

Furthermore, the erosion between the slope rills does not change continuously, because each 
grid element is connected to others in the x or y direction in the modelling process. The slope 
gully is susceptible to major pattern changes and substantial shifts as the slope gully is eroded, 
transported, and deposited, where the sediment load and water to discharge may change. The 
gradient in the y direction on each grid is steeper than the gradient in the x direction. Therefore, 
the soil is more easily eroded in the vertical direction, that is, in the y direction of the grid. 

As shown in Figure 3, the sediment yield and rill erosion width show an overall upward trend 
over time. The sediment yield and rill erosion width have very similar trends over time, and the 
time 800 s is the turning point when the rill erosion width begins to be greater than the sediment 
yield rate. The reason may be that the source of sand production on the slope has changed. The 
regression analysis results show that the cumulative sediment yield and rill erosion width have an 
extremely linear positive relationship: 

Y - aX — b, (12) 
where Y is the cumulative sediment yield (g); X is the rill erosion width (mm); and a and b are the 
undetermined coefficients, with the values of 0.8743 and 2.4012, respectively. In this study, a 
prototype rill is constructed along the flow direction in the middle of the slope to guide the 
development of the rill structure by assuming its width 1s 2.4012 mm. Therefore, in the case of 
different slope lengths and grades, modelling is carried out according to the actual situation. 
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Fig. 3 Relationship between sediment yield (Y) and rill erosion width (X) at different times 


As shown in Figure 4, the runoff velocity over time shows a downward trend; at test times of 
0—100, 100—140, and 320—510 s, the runoff velocity exhibits four sharp decreases. Two of the 
rapid decreases in the runoff velocity correspond to the two rapid increases in the narrow groove 
width. However, in the two phases at the beginning and end of the experiment, the changes in the 
flow rate of the concentrated water flow and the width of the rill are not consistent. In the first 
100 s of the test experiment, the runoff velocity 1s reduced by 40.096, while the width of the rill 
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remains unchanged. The main reason is that the shearing force of concentrated water on the ditch 
toe greatly increases the runoff width and reduces the runoff depth, so that the runoff velocity is 
reduced and the ditch foot wash-out height decreases with the increasing wash-out depth. Over 
320 to 360 s, the width of the narrow groove increases by 17.096, while the flow rate of 
concentrated water flow decreases by only 9.5%. Under the condition that the erosion rate of the 
ditch foot wash-out caused by the runoff shearing force is essentially unchanged, more ditch wall 
expansion and invasion occur. 
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Fig. 4 Variations of rill erosion velocity and rill erosion width with the passage of time 


The regression analysis results shown in Table 3 reflect that there is a significant negative 
correlation between runoff velocity and rill erosion width, which can be expressed by an 
exponential function: 


Y-104.12e 99? (R?=0.9591, P<0.005), (13) 


where Y is the runoff velocity (cm/s) and X is the width of the narrow groove (cm). 

The applicability and accuracy of the GAST (GPU Accelerated Surface Water Flow and 
Transport Model) in the study area can be determined by R^ and Eng indicators for evaluation. R^ 
is calculated by a linear regression between the measured value and the simulated value to 
represent the simulation consistency level of the change trend between runoff and measured flow. 
Eys represents the similarity between the simulated value and the measured values at the 1:1 level. 
For example, Exs-0.8 means that 80% of the simulated values are similar to the measured values. 
The calculation formula of Ens is as follows: 


Eys = 1- £—— ——;. (14) 


In the formula, M is the actual measured value, O is the model simulated value, M is the 


average of the actual measured values, O is the average of the simulated values, and n is the 
number of observation data points. 

As shown in Figure 5, the high-efficiency and high-resolution numerical model based on 
GPU-based acceleration technology that couples hydrodynamics and morphology established in 
this study has been indicated to be reliable and has achieved good results in simulation prediction. 
The calculation results show that the R? and Eys values are 0.97 and 0.98, respectively. Therefore, 
the erosion prediction model established in this study shows a strong and broad general type, 
which reflects that the model can well simulate the spatial distribution characteristics of 
inter-ditch erosion. 
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Table3 Erosion rill velocity and width at different times 


Time (s) Erosion rill velocity (cm/s) Erosion rill width (mm) 
1.36 75.20 1.74 
30.90 68.10 1:79 
58.10 61.50 4.3] 
93.00 53.40 10.30 
108.00 43.50 15.90 
145.00 38.60 18.40 
187.00 36.90 20.30 
236.00 34.70 21.20 
284.00 32.90 23.10 
339.00 28.30 26.00 
395.00 23.90 34.50 
470.00 17.80 40.70 
521.00 15.30 43.60 
584.00 12.70 45.90 
660.00 10.20 49.30 
742.00 9.02 53.90 
795.00 9.04 60.30 
858.00 8.20 65.00 
921.00 6.99 68.40 
966.00 6.89 72.20 
1000.00 7.02 76.90 


Simulated value 


0 5 10 15 20 25 30 35 40 45 50 
Measured value 
Fig. 5 Correlation between the simulated and measured rill erosion values 


4.2 Numerical simulation example of the erosion process in a small watershed 


The second case chosen to test this model is an experiment of Coquet River Basin that involves 
watershed erosion (Hou et al., 2020). The associated reaches are chosen because they have 
considerable historic instability in terms of channel patterns and morphology, and rising river 
water promotes erosion under rainfall conditions. As shown in Figure 6, the erosion depth 
changes and erosion cumulative amount of the basin over 1—40 h is simulated. The elevation 
difference of erosion is on the order of 0.5-3.0 m. The average erosion depth overall increases 
over time, especially at £—20—30 h (Fig. 6e-f). 
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Fig. 6 Numerical simulation of watershed erosion process of the Coquet River Basin at different times. (a), 1 h; 
(b), 5 h; (c), 10 h; (d), 15 h; (e), 20 h; (f), 30 h; (g), 40 h. 


For significant erosion sites and rainfall events, the key simulation results focus on determining 
model capabilities and erosion patterns. The results of these simulations appear as hourly outputs 
(data not shown). Selected important time frames are used to represent the changes in the spatial 
location of erosion over the entire basin. These results are used to illustrate the morphological 
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changes of watershed erosion, including sediment erosion and important river sections of 
deposition, to provide a concept of the relationship between morphological changes and eroded 
energy (Fig. 7). Due to river erosion and sediment exchange and migration (Ham, 2005), energy 
seems to concentrate in the riverbed (Fuller et al., 2003). In addition to the large amount of 
sediment migration that leads to river channel deposition, the river channel migrates again during 
the relaxation period due to erosion and then redeposits in the lower reaches of the river. The 
large amount of siltation that occurs in the river section of the river basin is also attributed to the 
increase in the inflow of sediment at the higher flow rate coupled with the input of river bank 
erosion, considering that the low slope of the study river section is a key factor that affects the 
river's response to major floods. This restricts the transport and deposition of coarse bedrock in 
the upper part (Danks, 1998). Therefore, a low slope leads to decreases in sediment transport 
capacity and water flow power. Accordingly, the high rate of bank retreat is related to the 
excessive concentration of flow energy on the bank, leading to greater lateral erosion, especially 
on the outside of the bend (Newson, 1980; Hickin, 1984). 


Erosion energy 


Fig. 7 Watershed erosion energy exchange. "+" and "—" indicate positive and negative work done, respectively. 
Positive work corresponds to erosion and negative work to deposition. 


As shown in Figure 8, a total of nine cross sections of the Coquet River Basin are considered 
for the measured old riverbed and simulated new eroded riverbed analysis. By analysing the 
conditions of each section, we can better understand the changes of the topography of the river 
section that are most affected by different erosion intensities and key deposition processes, as 
well as their influencing factors. Based on the results of a series of cross sections, the more stable 
subranges are more resistant to morphological changes and sediment flux, and these subranges 
have not yet undergone fundamental landform changes in response to high-intensity erosion. In 
terms of erosion sediment mobility, there are strong responses to high-intensity erosion or 
frequent low-amplitude events. 

As shown in Figure 8a, when comparing the old riverbed with the new eroded riverbed, at the 
distance of 20.00 m, the new eroded riverbed is lower than the old riverbed, with a height 
difference of approximately 0.45 m. This reduction is caused by processes such as strong erosion 
of the watershed and the "over-transport" of sediments. These processes change the structure of 
the upper reaches of the basin and transfer a large amount of sediment to the middle and lower 
reaches of the basin. Similarly, as shown in Figure 8b-i, the erosion change range of the old 
riverbed and the new eroded riverbed is 0.20—1.00 m. Continuous erosion processes and the 
increase in the width of the lateral soil erosion lead to the further migration of Cross section 9 
(Fig. 81), and bank retreat can potentially be driven by frequent bank failures on the left-hand 
bank. This process, also known as cantilever failure, is the progressive removal of riparian 
material by flow erosion of non-viscous riparian material, which is composed of non-viscous 
sand. 

The quantitative analysis results about the erosion and deposition processes of riverbed at the 
nine cross sections of the Coquet River Basin are shown in Figure 9. In the most severely eroded 
area, the erosion depth is approximately 13.00 m. In general, the degree of erosion for each cross 
section is greater than that of siltation. The depth of erosion change value ranges from —0.86 to 
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Fig. 8 Comparisons of elevations for the measured old riverbed and simulated new eroded riverbed at nine cross 
sections of the Coquet River Basin (a1) 


—2.79 m, and the depth of deposition change value is in the range of 0.38—1.02 m. Normally, the 
erosion depth of most areas is approximately 1.00 m, mainly in flood-prone regions. We compare 
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erosion and deposition by calculating the accumulation of sediment in this basin, and find that the 
overall erosion is net erosion. This illustrates that erosion process dominates the basin in terms of 
the erosion geomorphological bed and nine cross section profiles, in turn leading to an increase in 
the width of the basin. Therefore, the process can be regarded as watershed erosion. 

For this reason, these changes can be attributed to erosion and deposition, which result in 
alterations in the river channel course, bed elevation, and replacement due to sediment 
accumulation and changes in the riverbed gradient. In addition to the lateral development of the 
river course, the bends of the upstream river are eroded, increasing the depth of that particular 
section of the riverbed. The erosion depth at the Cross section 3 is 2.40 m. The erosion depth at 
the Cross section 4 is 1.29 m, and the erosion depth at the Cross section 6 is 0.63 m. As a result, 
the aforementioned factors clearly affect the active channel, particularly by altering the width and 
gradient of the channel, in addition to the consequential bank erosion. Thus, the eroded channels 
can be eroded into new channels during the aforementioned reciprocating erosion process and 
change their width and slope, thus enhancing the erosion process. 

Erosion and deposition of the riverbank characterize the basin, as determined by calculating the 
changes in the balance of erosion and deposition over the period for cross sections 1—9 (Fig. 9) in 
the basin. The net sediment gains and losses during this period indicate an excess of sediment 
erosion and the export of material downstream from the reach, which can also occur in stable 
rivers. This process is also linked to sediment distribution. Generally, the increase in sediment 
flux 1s associated with higher rainfall magnitude in conjunction with the sediment input of 
erosional bank processes, which lead to changes in watershed geometry. Estimating the riverbed 
and riverbank material erosion along the basin channel may provide an indication of the dominant 
processes that influence morphologic change and the development of the meandering basin. 
Consequently, if the eroded sediment is not removed or transported, then a change can be 
expected in the shape and form of a cross section related to the active basin width. 

To demonstrate the performance of the GPU code calculation efficiency, we use a model of 
CPU (CPU with four cores (Intel Xeon E5-2609/CPU)) calculation to simulate the same test case. 
The runtime uses different grid resolutions on different hardware devices. As the grid becomes 
finer, the performance of the GPU code in terms of computational efficiency improves. Therefore, 
in soil erosion events, GPU acceleration can achieve large-scale simulation calculations and play 
a great role in the evaluation of water and soil loss processes. 


5 Conclusions 


This study develops a new soil erosion model that couples 2D hydrodynamics and the G-A 
rainfall infiltration model. The model also applies GPU-based acceleration technology to increase 
the calculation speed. The model tested shows that its simulation values and measured values 
indicate a good consistency, and the slope erosion patterns and spatial distribution characteristics 
at different times are highly consistent with the final ravine shape. It is possible to make good 
predictions of slope erosion using this developed model. A linear function can be used to predict 
the amount of sediment transport and the width of rill erosion in slope soil erosion. The rate of 
erosion depends mainly on the resistance of the soil to corrosion. There are two turning points in 
the erosion rate and rill erosion width. When the velocity is sharply reduced, the erosion width 
gradually increases, and when the two reach an equilibrium point, the erosion width gradually 
decreases as the velocity increases. Notably, the width of the erosion rill depends mainly on the 
formation and development of the groove opening from small to large. The effect of the shearing 
force, although small, cannot be ignored. In large-scale basins, bank incisions collapse from small 
to large, resulting in a large area of riverbank to be eroded and washed away. 


Acknowledgements 


This research was funded by the National Natural Science Foundation of China (52079106, 52009104, 51609199) 
and the National Key Research and Development Program of China (2016YFC0402704). We also thank the 
anonymous reviewers and editors for their constructive comments for this manuscript. 


KANG Yongde et al.: Two-dimensional hydrodynamic robust numerical model of soil erosion... 


NS 

O 

— 

[^] 

= 

i] 

- 

un 

E 
N 
fæl 
9 
5 
o 
ua 
un 
un 
E 
e 
nn 
= 

t a + 

= TT 


Distance (m) 


oss section 1 


NSSSSSSSSNNS 
NQQQ QW °” Fe 
KX gg). 


(a) Cr 


SS Si oe a 
et ts ia 


(ur) ydəp uorisodoqqgdop uorsorq 


SS 


Distance (m) 


SN 


SESS ETT S D 
(ur) ydəp uonisodoqq/ydop uorsoxq 


NSW WC, 
MMW GG 


Wry 


Distance (m) 


(c) Cross iar 3 


XM ME M. mee eir 
S S SoSo 


(w) ydəp uonrsodoqq/qgidop uorsorq 


1o 
Wu 


(f) Cross section 6 
Distance (m) 


0.2 


MW WW) 


Ga 


Distance (m) 


A 
WW 


IW GG 
.W,,. WW 


(e) uu section 5 


-].0 


a + 
TH 


(w) dep uonrsodo(q/qidop uorsorq 


SS 


NSW 


= 


(h) Cross section 8 


(g) Cross section 7 
- 
- 


(w) undop uonisodo(qq/ugdop uorsorq 


Distance (m) 


23 45 67 8 sane 
Distance (m) 


RSs SSS 

SSSSSSNS 

SSSSSSSSSSSSSSNSNSSNNNNNNNN 

SX )BW 
SSS AN 


Rew 


S ~ 

>on g 

ENSNSN- A E 

NS 
AA n d 
SSSSSSSSSSSSSSSSSSSSSNWNNSNSNS z 
SH Way ES 
SW) + 
SSN A 
RAAT A 
SSX Ss 

Os SSSSSSSSSSSSSSSSNNS 

e [SSSSSSSSSSSSSSS 

o SSSSSSSSNS 

HG SSSsg 

oO N| 

D ISSN 

na SS 

n SW ws 

a SS SSS 

o RQ MTT 

H SSS NN 

= N 
RN 
ANN 

NO SOOO Q xo 00 O C xo CO 

vier edes SIN D e 


(ur) ydəp uonisodoqq/gdop uorsorq 
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